Method, computer program, and system for intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale

ABSTRACT

A method and system for intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale including receiving an input signal, determining a baseline segment and a monotonic residual segment with strictly negative minimum and strictly positive maximum between two successive extrema of the input signal, and producing a baseline output signal and a residual output signal. The method and system also includes determining at least one instantaneous frequency estimate from a proper rotation signal, determining a zero-crossing and a local extremum of the proper rotation signal, and applying interpolation thereto to determine an instantaneous frequency estimate thereof. The method and system further includes determining at least one instantaneous frequency estimate from a proper rotation signal, extracting an amplitude-normalized half wave therefrom and applying an arcsine function to the amplitude-normalized half wave to determine an instantaneous frequency estimate of the proper rotation signal.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation application of U.S. application Ser.No. 10/684,189 filed Oct. 10, 2003, now U.S. Pat. No. 7,054,792, issuedMay 30, 2006; which claims priority of Provisional Patent ApplicationNo. 60/418,141 filed Oct. 11, 2002.

COMPUTER PROGRAM LISTING APPENDIX

An appendix containing five computer files on compact disk is includedin this application. The files, which are formated for an IBM PC/MSWindows-compatible computer system, are as follows: contents.txt 864bytes; itd.m 964 bytes; itd_step.m 1,528 bytes; itd_sift.m 1,477 bytes;pangle.m 2,881 bytes; Itd.cpp 3,619 bytes; Itd.h 1,742 bytes;ItdStep.cpp 12,300 bytes; ItdStep.h 6,180 bytes; itd_rec.cpp 3,384bytes; FhsBuffers.h 4,275 bytes;These files were loaded on the non-rewritable CD-R disc on Oct. 7, 2003and are incorporated herein by reference.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to a system for analyzing signals anddata and, more specifically but without limitation, to a system foranalyzing non-linear and/or non-stationary signals and data.

2. Description of Related Art

Data analysis is an essential part of pure and applied research in manydisciplines. In many practical applications, the raw data havecharacteristics that change over time, commonly referred to as beingnon-stationary and, additionally, very often result from non-linearprocesses. Unfortunately, the majority of existing methods for analyzingdata are designed to treat stationary, linear data.

Another common and serious problem of data analysis is the existence ofnoise and/or non-stationary trend. Common practice to deal with theseproblems has involved application of band pass filters to the data.However, these filters are Fourier based and, as such, typically resultin the introduction of spurious harmonics in non-stationary data.Therefore, Fourier-based filters have limited utility and practicalvalue for use with non-stationary and/or non-linear data. In addition,the (low frequency) signal trend can carry significant and usefulinformation about the process being analyzed, and thus should not besimply filtered out. Prior art methods for processing non-stationarydata include Fourier analysis, wavelet analysis, the Wigner-Villedistribution, the evolutionary spectrum, the empirical orthogonalfunction expansion, and the empirical mode decomposition. These priorart methods can be briefly described as follows.

Fourier Analysis.

Traditional methods of time-frequency-energy analysis are based onFourier transforms and are designed for stationary and linear processes.The application of these methods to analysis of non-stationary,non-linear data can give misleading results. For example, the Fourierspectrum defines uniform harmonic components globally. Therefore,Fourier analysis needs many additional harmonic components in order tosimulate non-stationary data, which are not uniform globally. As aresult, Fourier analysis tends to spread signal energy over a widefrequency range. As is well known by those having skill in the art, thefaster the change in the time domain, the wider the frequency range.Unfortunately, many of the components, that are added in order tosimulate the non-stationary nature of the data in the time domain,divert energy to a much wider frequency domain. For example, a singleimpulse, a signal whose deviation from constancy occurs at a singlemoment in time, requires infinitely many frequencies with identicalpower to represent it. Constrained by energy conservation, the spuriousharmonics that are added and the wide frequency spectrum required tosimulate the non-linearity cannot faithfully represent the true energydensity in the resulting frequency space.

Further, Fourier spectral analysis utilizes linear superposition oftrigonometric functions. Such analysis needs additional harmoniccomponents in order to simulate the effects of non-linearity, such asdeformed wave profiles. Whenever the form of the data deviates from apure sine or cosine function, the Fourier spectrum will containharmonics. As explained above, both non-stationarity and non-linearitycan induce spurious harmonic components that cause energy spreading. Theresulting consequence is a misleading or incorrect time-frequencydistribution of signal energy for non-linear and/or non-stationary data.

Many data analysis methods have been developed based on Fouriertransforms. The spectrogram is the most basic and common method, whichis a finite-time window Fourier spectral analysis that is repeated inmoving-time windows. By successively sliding the window along the timeaxis, a time-frequency-energy distribution is obtained. Since such adistribution relies on the traditional Fourier spectral analysis, themethod assumes the data to be piecewise stationary. This assumption isnot valid for most non-stationary data. Even if the data were piecewisestationary, it is highly unlikely in most cases that the window sizeadopted would coincide with the stationary time scale. Furthermore,there are practical difficulties in applying the method. In order tolocalize an event in time with good temporal precision, the window widthmust be narrow. Unfortunately, the frequency resolution worsens aswindow size decreases. Although the conflict in these requirements canbe mitigated by different techniques, it renders the applicability ofFourier analysis to non-linear, non-stationary data of limited use.

Wavelet Analysis.

Wavelet analysis, which has become extremely popular during the pastdecade, is an attempt to overcome the problems of windowed Fourieranalysis by utilization of a basis of functions to represent a signalthat contains elements having different time scales. This approachallows wavelet analysis to detect changes that occur rapidly, i.e.,those on a small time scale, with good temporal resolution but poorfrequency resolution, or slowly, i.e., those on a large time scale, withgood frequency resolution but poor temporal resolution. Morespecifically, the wavelet analysis approach is essentially anadjustable-window Fourier spectral analysis. Wavelet analysis is usefulfor analyzing data with gradual frequency changes. Primary applicationsof wavelet analysis have been in areas of edge detection and audio andimage signal compression. Limited applications also include analysis oftime-frequency distribution of energy in time series and oftwo-dimensional images such as fingerprints. Unfortunately, the finitelength of the basic wavelet function results in energy leakage acrossdifferent levels of resolution in a multi-resolution analysis, whichcauses quantitative applications of the time-frequency-energydistribution to be more difficult.

Sometimes, the interpretation of the wavelet can also becounterintuitive. For example, the more temporally localized the basicwavelet, the higher the frequency range will be. Therefore, to define achange occurring locally, the analytic result may occur in a highfrequency range. In other words, if a local event occurs only in a lowfrequency range, the effects of that local event may only appear in ahigh frequency range. In many applications, the interpretation of such aresult would be difficult if not impossible.

Another difficulty with wavelet analysis is that it is not adaptive.Once the basic wavelet is selected, the basis for the analysis iscompletely determined to the extent obtainable from the selected basicwavelet, and all information of the input signals is represented interms of that basis. Although the basis can be specially selected for anindividual application, the information obtained depends heavily on theproperties inherent to that basis rather than solely on the intrinsicproperties of the signals being studied. Malvar wavelets, WaveletPackets, and Matching Pursuit methods have been developed to overcomesome of these limitations to more accurately represent a signal havingdynamics that vary with time and that include both stationary andnon-stationary characteristics. Unfortunately, these developments inwavelet analysis continue to suffer from the representation of thesignal information in terms of a pre-selected basis of functions thatoften has little or nothing to do with the dynamics and othercharacteristics of the input signals.

The Wigner-Ville Distribution.

The Wigner-Ville distribution, sometimes referred to as the Heisenbergwavelet, is the Fourier transform of the central covariance function.The difficulty with this method is the severe cross terms indicated bythe existence of negative power for some frequency ranges. TheWigner-Ville distribution has been used to define wave packets thatreduce a complicated data set into a finite number of simple components.Although this approach can be applied to a wide variety of problems,applications to nonstationary or nonlinear data are extremelycomplicated. Further, such applications again suffer from the samelimitation of the other prior art methods described above in that thebases for representation of information are not derived from the dataitself.

Evolutionary Spectrum.

The Evolutionary Spectrum approach is used to extend the classic Fourierspectral analysis to a more generalized basis, namely from sine orcosine functions to a family of orthogonal functions indexed by time anddefined for all real frequencies. Thus, the original signal can beexpanded into a family of amplitude modulated trigonometric functions.The problem with this approach, which severely limits its applicability,is due to the lack of means for adequately defining the basis. Inprinciple, the basis has to be defined a posteriori in order for thismethod to work. To date, no systematic method of constructing such abasis exists. Therefore, it is impossible to construct an evolutionaryspectrum from a given set of data. As a result, applications of theevolutionary spectrum method have changed the approach from dataanalysis to data simulation. Thus, application of the evolutionaryspectrum approach involves assumptions causing the input signal to bereconstituted based on an assumed spectrum. Although there may be somegeneral resemblance of a simulated input signal to the correspondingreal data, it is not the data that generated the spectrum. Consequently,evolutionary spectrum analysis has very little useful application, ifany.

The Empirical Orthogonal Function Expansion.

Empirical orthogonal function expansion (“EOF”), also known as“principal component analysis” or “singular value decomposition,”provides a means for representing any function of state and time as aweighted sum of empirical eigenfunctions that form an orthonormal basis.The weights are allowed to vary with time. EOF differs from the othermethods described hereinabove in that the expansion basis is derivedfrom the data. The critical flaw of EOF is that it only gives a variancedistribution of the modes defined by the basis functions, and thisdistribution by itself gives no indication of time scales or frequencycontent of the signal. In addition, any single component of thenon-unique decomposition, even if the basis is orthogonal, does notusually provide any physical meaning. The Singular Spectral Analysis(“SSA”) method, which is a variation of EOF, is simply the Fouriertransform of EOF. Unfortunately, since EOF components from a non-linearand non-stationary data set are not linear and stationary, use ofFourier spectral analysis generated by SSA is flawed. Consequently, SSAis not a genuine improvement of EOF. Although adaptive in nature, theEOF and SSA methods have limited applicability.

The Empirical Mode Decomposition.

The empirical mode decomposition (“EMD”) method, involves two majorsteps. The first step is the application of an algorithm for decomposingphysical signals, including those that may be non-linear ornon-stationary, into a collection of Intrinsic Mode Functions (“IMFs”),which are supposedly indicative of intrinsic oscillatory modes in thephysical signals. More specifically, the cornerstone of the EMD methodis the extraction of a signal baseline from a physical signal whereinthe baseline is computed as the mean value of the upper and lowerenvelopes of the physical signal. The upper envelope is defined by cubicsplines connecting the local maxima of the physical signal, and thelower envelope is defined by cubic splines connecting the local minimaof the physical signal. The signal baseline is then extracted, orsubtracted, from the original signal to obtain an IMF having the firstand highest frequency present in the signal.

A goal of the first step is to obtain a well-behaved IMF prior toperforming the second step: applying a Hilbert Transform to the IMF.“Well-behaved” means that the IMF should be a “proper rotation,” i.e.,all local maxima are strictly positive and all local minima are strictlynegative. This does not necessarily happen with one step of EMD, andthus a laborious, iterative “sifting” process is applied to the signalbaseline and is terminated when a set of “stopping criteria” aresatisfied, such as when the resulting IMF either becomes a properrotation or when some other criteria are reached (e.g., a predeterminednumber of iterations exhausted) without obtaining a proper rotation. Thestopping criterion is based on a combination of (i) limiting the amountof computational energy expended, and (ii) having the constructed IMFclosely approximate the desired property. When the first IMF functionhas been obtained, it is subtracted from the signal and the process isrepeated on the resulting lower frequency signal. This process isrepeated again and again until the decomposition is completed in thewindow of signal being analyzed.

The sifting process has two goals: (i) to separate out high frequency,small amplitude waves that are “riding” atop, or superimposed on, largeramplitude, lower frequency waves, and (ii) to smooth out unevenamplitudes in the IMF being extracted. Unfortunately, these goals areoften conflicting for non-stationary signals wherein riding waves may beisolated and/or are highly variable in amplitude. As a result, thesifting process must be applied cautiously as it can potentiallyobliterate the physically meaningful amplitude fluctuations of theoriginal signal.

As noted, once the IMFs have been obtained, the second step of the EMDmethod is to apply the Hilbert Transform to the IMFs which, provided theIMFs are well-behaved, results in quantified instantaneous frequencymeasurements for each component as a function of time.

However, the EMD method suffers from a number of shortcomings that leadto inaccuracies in depiction of signal dynamics and misleading resultsin subsequent signal analysis, as follows:

(a) The construction of the IMF baseline as the mean of the cubic splineenvelopes of the signal suffers from several limitations, including thetime scale being defined only by the local extrema of the originalsignal, ignoring the locations of the rest of the critical points, suchas inflection points and zero crossings, which are not preserved by thesifting process.

(b) The EMD transformation is window-based; the sifting procedure andother processing requires an entire window of data to be repeatedlyprocessed. The resulting information for the entire window is notavailable until the window's processing has been completed. On average,this causes a delay in obtaining the resulting information of at leasthalf the window length plus the average computational time for a windowof data. As a result, an insurmountable problem is created for the EMDin real-time/online applications. In order to obtain the informationquickly, the window length should be short. However, short windows haveless accurate results due to boundary or edge effects, and are incapableof resolving frequency information on a time scale longer than thewindow length itself.

(c) The EMD method is computationally expensive, and also subject touncertain computational expenses. The EMD method requires repeatedsifting of components in order to obtain well-behaved IMFs or untilstopping criteria are satisfied. Such a procedure may require numerousiterations and may not occur in finite time. Thus, the method often willnot result in IMFs with the desired proper rotation property, even whensifting numerous times. Also, the IMFs are generally dependent on theparameters of the algorithm that define the stopping criteria forsifting.

(d) Overshoots and undershoots of the interpolating cubic splinesgenerate spurious extrema, and shift or exaggerate the existing ones.Therefore, the envelope-mean method of the extraction of the baselinedoes not work well in certain cases, such as when there are gaps in thedata or data are unevenly sampled, for example. Although this problemcan be mitigated by using more sophisticated spline methods, such as thetaut spline, such improvements are marginal. Moreover, splines are notnecessarily well suited to approximate long timescale trends in realdata.

(e) The EMD method does not accurately preserve precise temporalinformation in the input signal. The critical points of the IMFs, suchas the extrema, inflection points, etc., are not the same as those ofthe original signal. Also, the EMD method, being exclusively determinedby extrema, is deficient in its ability to extract weak signals embeddedin stronger signals.

(f) Since the cubic splines can have wide swings at the ends, theenvelope-mean method of the EMD method is particularly unsuitable forreal-time applications or for applications utilizing a narrow window.The end effects also propagate to the interior and significantly corruptthe data, as the construction of IMFs progresses, as can be seen inFIGS. 1-2. FIG. 1, left panel, illustrates the application of the EMD toa test signal (top-most signal) to produce IMF components (displayedbelow test signal). This panel illustrates end effects, spline-relatedinstabilities (most noticeable in bottom components), and inability toextract the readily apparent signal baseline. The intrinsic timescaledecomposition, sometimes referred to herein as ITD, (shown in the rightpanel) separates the same test signal (top-most signal) into stablecomponents (displayed below test signal), demonstrating the fact thatthat ITD has no end effect propagation beyond the first two extrema ateach level and allows correct identification of the trend (dashed line).FIG. 2, top panel, shows a brain wave input signal (electrocorticogram;‘ECoG’) containing an epileptic seizure used to illustrate decompositiondifferences between EMD and ITD. EMD-generated (lower left panel) andITD-generated (lower right panel) decompositions of the cumulative sumof the raw signal show that ITD, unlike EMD, does not generateextraneous components and correctly reveals large timescale variationsof the signal (DC trend).

(g) The application of the Hilbert Transform to track instantaneousfrequency is only appropriate when frequency is varying slowly withrespect to amplitude fluctuations which, unfortunately is a conditionnot satisfied by many non-stationary signals. Moreover, properrotations, which are not guaranteed by the EMD, are necessary for theexistence of a meaningful, instantaneous frequency that is restrictiveand local.

(h) Primarily due to the sifting procedure, the EMD method causes (i) asmearing of signal energy information across different decompositionlevels, and (ii) an intra-level smoothing of energy and frequencyinformation that may not reflect the characteristics of the signal beingdecomposed or analyzed. The potential negative effects of over-siftinginclude the obliteration of physically meaningful amplitude fluctuationsin the original signal. Also, the inter-level smearing of energy andlimitation of decomposition in a window of data create inaccuracies inthe EMD's representation of the underlying signal trend, especially ifthe signal trend is longer than a single window.

(i) The EMD method cannot guarantee that the IMF components will be“proper rotations,” even with the sifting reiterations. Often, inpractice, they are not.

(j) The EMD method is not well behaved if the data is unevenly sampledor if it is discontinuous, and, therefore, may not preserve intact phaseand morphology characteristics of the signal.

What is needed is a system for reliably analyzing non-linear and/ornon-stationary signals and data capable of decomposing them and/orextracting precise time-frequency-energy distribution information.

SUMMARY OF THE INVENTION

In a method for decomposing an input signal into a baseline signal and aresidual signal, the steps including receiving an input signal into aprocessor, determining a monotonic residual segment with strictlynegative minimum and strictly positive maximum and a baseline segmentwherein said segments are defined on a time interval comprising theinterval between two successive extrema of the input signal and whereinthe input signal on that interval is the sum of the baseline andresidual segments, and producing a baseline output signal and a residualoutput signal wherein the baseline signal is obtained from the baselinesegment and the residual signal is obtained from the residual segment asdetermined, such that the sum of the baseline and residual signals isequal to the input signal thereby forming a decomposition of the inputsignal.

In a method for signal decomposition, the steps including receiving aninput signal into a processor, using the input signal to construct abaseline signal component such that subtraction of the baseline signalcomponent from the input signal always produces a proper rotationsignal.

In a method for determining at least one instantaneous frequencyestimate from a proper rotation signal, the steps including inputing aproper rotation signal to a processor, determining a zero-crossing and alocal extremum of the proper rotation signal, and applying linearinterpolation to the proper rotation signal between the zero-crossingand the local extremum to determine an instantaneous frequency estimateof the proper rotation signal.

In a method for determining at least one instantaneous frequencyestimate from a proper rotation signal, the steps including inputing aproper rotation signal to a processor, extracting anamplitude-normalized half wave from the proper rotation signal, andapplying an arcsine function to the amplitude-normalized half wave todetermine an instantaneous frequency estimate of the proper rotationsignal.

PRINCIPAL OBJECTS AND ADVANTAGES OF THE INVENTION

The principal objects and advantages of the present invention include:providing a system and method for analyzing non-stationary and/ornon-linear signals; providing such a system and method that can operatein real-time due to its computational simplicity, recursive nature,elimination of need for sifting, and absence of significant end effects;providing such a system and method that adapts to any timescale and usescomplete signal information, including all critical points such asinflection points and zero crossings; providing such a system and methodthat can extract weak signals embedded in stronger signals; providingsuch a system and method that provides precise time-frequency-energylocalization, either simultaneously via the Hilbert transform or novelmethods for instantaneous phase and instantaneous frequency computationnot requiring use of this transform, on other timescales, such as theinter-extrema timescale or the inter-critical point timescale forexample; providing such a system and method having the ability tointerpolate adjacent critical points using the signal itself; providingsuch a system that is well-behaved even if the data is unevenly sampledor is discontinuous; providing such a system that may be applied toanalysis of multi-dimensional signals or data such as images orsurfaces; providing such a system and method that has the ability topreserve intact phase and morphology characteristics of the signal;providing such a system and method that can extract proper rotationcomponents in one step; providing such a system and method wherein allextracted components are guaranteed to be “proper rotations” with allstrictly positive local maxima and all strictly negative local minima;providing such a system and method wherein the decomposition iscompletely determined by the input signal; providing such a system andmethod that is fully adaptive and local in time; providing such a systemand method wherein each step consists only of comparisons to determineextrema followed by a piece-wise linear transformation of buffered databetween two successive extrema to produce the desired signal components;providing such a system and method that preserves temporal informationin the signal; providing such a system and method wherein the criticalpoints of all extracted components, such as the points in time at whichlocal maxima, local minima, inflection points, etc., occur, coincideprecisely with critical points of the original signal; providing such asystem and method that allows waveform feature-based discrimination tobe used in combination with single wave analysis and classification toproduce powerful and flexible new signal filters for use indecomposition, detection and compression; providing such a system andmethod that substantially eliminate boundary or windowing effects;providing such a system and method for determining at least oneinstantaneous frequency estimate from a proper rotation signal usinglinear interpolation of the proper rotation signal between azero-crossing and a local extremum; providing such a system and methodfor determining at least one instantaneous frequency estimate from aproper rotation signal using an arcsine function applied to anamplitude-normalized half wave of the proper rotation signal; andgenerally providing such a system and method that is effective inoperation, reliable in performance, capable of long operating life, andparticularly well adapted for the proposed usages thereof.

Other objects and advantages of the present invention will becomeapparent from the following description taken in conjunction with theaccompanying drawings, which constitute a part of this disclosure andwherein are set forth exemplary embodiments of the present invention toillustrate various objects and features thereof.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A-1B is a comparison of the prior art method of EMD with the ITDas applied to a test signal.

FIGS. 2B-2C is a comparison of the prior art method of EMD with the ITDas applied to a brain wave signal as shown in FIG. 2A.

FIG. 3 is an illustration of ITD's extraction of the baseline from aninput signal.

FIGS. 4A-4D are illustrations of the steps of the online ITD method.

FIGS. 5B-5D are comparisons between the prior art methods of Fourieranalysis and wavelet analysis and the ITD (FIGS. 5E-5J) in determiningtime-frequency-energy (TFE) distributions from a sample signal as shownin FIG. 5A.

FIGS. 6A-6G are illustrations of the ITD-based method.

FIGS. 7A-7H are illustrations of the ability of the ITD-based filteringmethod, used to differentiate between two types of waves that havesignificantly overlapping spectral characteristics.

FIGS. 8A-8C are illustrations of ITD applied to multidimensional data.

FIG. 9A is a flow chart representation of the ITD method of the presentinvention.

FIG. 9B is a detailed flow chart representation of a portion of the ITDmethod of the present invention as indicated at 9B of FIG. 9A.

FIG. 9C is a detailed flow chart representation of a portion of the ITDmethod of the present invention as indicated at 9C of FIG. 9A.

FIG. 9D is a detailed flow chart representation of a portion of the ITDmethod of the present invention as indicated at 9D of FIG. 9A.

FIG. 9E is a detailed flow chart representation of a portion of the ITDmethod of the present invention as indicated at 9E of FIG. 9A.

FIG. 9F is a detailed flow chart representation of a portion of the ITDmethod of the present invention as indicated at 9F of FIG. 9A.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

As required, detailed embodiments of the present invention are disclosedherein; however, it is to be understood that the disclosed embodimentsare merely exemplary of the invention, which may be embodied in variousforms. Therefore, specific structural and functional details disclosedherein are not to be interpreted as limiting, but merely as a basis forthe claims and as a representative basis for teaching one skilled in theart to variously employ the present invention in virtually anyappropriately detailed structure.

Briefly stated, the present invention comprises a system for automateddecomposition and analysis of signals of arbitrary type, origin, ortimescale. It is able to accurately analyze complex signals that may be,for example, non-stationary, nonlinear, or both. The decompositionobtained by the system of the present invention is completely adaptiveto the timescale of the analyzed signal (i.e., the input to the ITD), asdetermined by critical points of the input, such as local extrema, forexample. This system can further analyze the signal components resultingfrom this decomposition to accurately quantify and localize varioussignal characteristics in time and frequency. Such signalcharacteristics include amplitude, wavespeed, phase, regularity,morphology, moments, energy, variance, skewness, and kurtosis, forexample. The system can also apply such quantification and localizationinformation as a new type of adaptive filtering, such as for noiseremoval, trend extraction, or decomposition of the input signal intodifferent components, which have various desired properties. The ITD maybe also used as a real-time measure of signal complexity or ofinformation content as a function of time (or state) as determined bythe number of extractable components and the energy content existent ineach of them. For example, if a signal sample in a particular timewindow is decomposed by the ITD into 5 proper rotation components and atrend, and another time window of a (same or different) signal isdecomposed by the ITD into 9 proper rotation components and a trend,this indicates that the latter signal had more riding waves and thus mayhave carried more “information” than the former.

A. The ITD Method

The present invention includes an improved data analysis system foranalyzing non-stationary or non-linear signals or data. This method isbased on an algorithm, sometimes referred to herein as the IntrinsicTimescale Decomposition (“ITD”) algorithm, which decomposes a signalinto a set of components having sequentially lower frequencycharacteristics, according to the signal's intrinsic timescale asdetermined by local extrema, maxima and minima, or other criticalpoints. These timescales are important in quantifying and analyzing anoscillating or fluctuating system. The separation process relies on a“signal-based,” self-extraction, interpolation procedure used toidentify and isolate individual waves (it can also separate parts ofwaves and/or sets of waves). This process is then iterated wherein eachstep separates out the smallest remaining timescale modes (i.e., highestfrequency component) from the next successively larger timescale (i.e.,lower frequency) baseline signals. The iteration steps are repeateduntil all modes have been extracted, with the last such modecorresponding to the largest timescale present in the signal, sometimesreferred to herein as the monotonic signal “trend.”

This decomposition process results in a set of signal components which,if recombined, would sum to the original signal and which possess anumber of important properties that facilitate further signal analysis.These components contain constituent monotonic segments, each havingfully preserved temporal information (including precise preservation oftemporal location of all critical points of the original signal), thatcan be further decomposed, individually analyzed, or reassembledaccording to “wave speed,” energy, morphology, probability ofoccurrence, temporal localization, or any other feature of the segmentthat may be quantified. These structural elements provide a basis for anexpansion of the raw data in terms of the data itself. Most importantly,this basis is adaptive to arbitrarily fast or slow changes in theamplitude and/or timescale of the signal, which makes it ideallysuitable for analysis of non-stationary data.

One skilled in the art will recognize that it may be desirable incertain applications for a system based on ITD, as an operator onfunctions, to be applied to a prior transformation of a given functionand possibly followed by subsequent transformation (e.g., the inverse ofthe prior transformation). For example, a system based on ITD can beapplied to a signal after first differentiating it one or more times,after which the results may be integrated the same number of times toproduce the desired signal components. This technique can be especiallyuseful in uncovering embedded high frequency signals that have a verysmall signal-to-noise ratio.

A Single Step of the Method of the Present Invention.

The system of the present invention for extracting the baseline can beapplied to unevenly sampled data and preserves the locations of allcritical points of the original signal, i.e., points at which somesignal derivative of a given order is zero. FIG. 3 illustrates theprocedure involved in a single step of the method.

The height of a single lobe of the signal can be defined as the lengthof a vertical line drawn from the k^(th) extremum of the originalsignal, x(t), to the straight line connecting the extrema at x(t_(k−1))and x(t_(k+1)), and is given by$h_{k} = {{x\left( t_{k} \right)} - \frac{{\left( {t_{k} - t_{k - 1}} \right){x\left( t_{k + 1} \right)}} + {\left( {t_{k + 1} - t_{k}} \right){x\left( t_{k - 1} \right)}}}{t_{k + 1} - t_{k - 1}}}$where t_(k) is the time of the k^(th) extremum of x(t). The k^(th)segment of the baseline is then defined byb(t)=c _(k) +d _(k) x(t) for t _(k) ≦t≦t _(k+1),where c_(k) and d_(k) are constants. By requiring the values of thebaseline at the points, t_(k), to be b(t_(k))=x(t_(k))−h_(k)/2, then thetotal baseline can be constructed by concatenating these segments. Thevalues of c_(k) and d_(k) can then be determined by equating the valuesof the baseline at the junction points, namelyc _(k) +d _(k) x(t _(k))=x(t _(k))−h _(k)/2andc _(k) +d _(k) x(t _(k+1))=x(t _(k+1))−h _(k+1)/2.

The first iteration extracts the component r(t)=x(t)−b(t), whichcontains the smallest timescale oscillations present in the originalsignal, leaving larger timescale modes behind in the remaining“baseline”, b(t). The following properties then hold:

-   -   1. Each point b(t_(k)) is either an extremum or inflection point        with zero first derivative.    -   2. If x(t_(k+1))−x(t_(k))=(h_(k+1)−h_(k))/2, the k^(th) segment        (for t_(k)≦t≦t_(k+1)) of the baseline is a constant, b(t)=c_(k).    -   3. The component, r(t), is a proper rotation. That is, each        maximum of r(t) has a strictly positive value, and each minimum        has a strictly negative value. In addition, each point r(t_(k))        is either an extremum or an inflection point with a zero first        derivative and the locations of all other critical points are        preserved.

The baseline function and the extracted r(t) component are sometimesreferred to herein as LF(t) and HF(t), respectively, to indicate theirrelative frequencies as “low” and “high,” respectively. FIG. 3illustrates the baseline signal obtained by this single-step of themethod for a given input signal (x(t), also shown).

Iteration to Produce Multiple Levels

After the original signal has been decomposed into a high frequencycomponent, HF(t), and a lower frequency baseline component, LF(t), thisbaseline component can be treated as the original signal and similarlydecomposed in another level of the decomposition. Repetitive extractionof each successive baseline as the decomposition level increasesconverts: (i) more and more extrema into inflection points while makingthe lobes wider, and (ii) more and more of the data into long, monotonicsegments. The remainder of such a repetitive baseline extraction will bethe “trend,” a monotonic segment with timescale equal to the length oforiginal signal being analyzed.

The present invention, sometimes referred to herein as the IntrinsicTimescale Decomposition System, or the ITD system, includes an algorithmfor decomposing an input signal into components with successively lowerfrequency characteristics, until finally one is left with either amonotonic signal trend if all modes have been extracted and no extremaof the resulting signal remains, or a signal of lowest relativefrequency if the user limits the process to a fixed number of componentswhich are reached before the monotonic trend has been obtained. Themethod overcomes all of the limitations of prior art methods andprovides a significant advancement in the state-of-the-art.

In an application of the ITD method in order to decompose and analyze aninput signal, the following steps are followed:

1. The input signal, x(t), is divided into two component signals: (i) ahigh frequency component signal, sometimes referred to herein as the HFcomponent, that, adaptively at each point in time, contains the highestfrequency information present in the signal, and (ii) a baselinecomponent, sometimes referred to herein as the LF component or lowercomponent, that contains all remaining, relatively lower frequencyinformation of the signal. Thus, this first decomposition can bedescribed by the following equation:x(t)=HF(t)+LF(t)

2. The same process that was used to decompose x(t) is next applied tothe LF(t) “baseline” component that resulted from the first level ofdecomposition. In other words, LF(t) is treated as if it were theoriginal signal. By so doing, the LF(t) component is further subdividedinto two more components: (i) a component comprising the highestfrequency information of the LF(t) component, sometimes referred toherein as the HLF(t) component, and (ii) another component that containsall remaining relatively lower frequency information, sometimes referredto herein as the L²F(t) component. Thus, this second decomposition isdescribed by the equation:LF(t)=HLF(t)+L ² F(t).

3. Next, the L²F(t) component is similarly decomposed intoHL²F(t)+L³F(t), and the process is iterated on the successive L^(n)F(t)signals that contain successively lower and lower frequency information.At each step, the highest frequency information remaining is extractedor separated out. The process terminates when either an increasing ordecreasing monotonic L^(n)F(t) segment has been obtained (sometimesreferred to herein as the signal “trend”), or when the decomposition hasbeen performed a desired number of times. The resulting decompositioncan be summarized by the following equations: $\begin{matrix}{{x(t)} = {{H\quad{F(t)}} + {L\quad{F(t)}}}} \\{= {{H\quad{F(t)}} + {{HLF}(t)} + {L^{2}{F(t)}}}} \\{= {{H\quad{F(t)}} + {{HLF}(t)} + {{HL}^{2}{F(t)}} + {L^{3}{F(t)}}}} \\{= {{H\quad{F(t)}} + {{HLF}(t)} + {{HL}^{2}{F(t)}} + {{HL}^{3}{F(t)}} + {L^{4}{F(t)}}}} \\{= {{H\quad{F(t)}} + {{HLF}(t)} + {{HL}^{2}{F(t)}} + {{HL}^{3}{F(t)}} + \ldots +}} \\{{{HL}^{n - 1}{F(t)}} + {L^{n}{F(t)}}}\end{matrix}$

As the decomposition is being performed, each new extrema in one of thecomponents generates a new monotonic segment in the next lower frequencylevel, with properties that can be easily quantified and analyzed. Asthe original signal is decomposed by the method into components, it issimultaneously broken down into an ensemble of individual monotonicsegments that are locked in time with the original extrema of thesignal, i.e., these segments span the time interval between localextrema of the original signal These individual segments can beclassified according to their own characteristics, allowing new anduseful filtered signals to be constructed by assembling those segments,lobes, waves or groups of waves that have certain specificcharacteristics or properties. This process allows the invention to beapplied to create an entirely new type of non-linear signal filter thatis able to differentiate precise time-frequency-energy informationsimultaneously with many other important wave-shape characteristics atthe single wave level in a computationally efficient manner.

Such decomposition, as provided by the ITD system of the presentinvention, has a number of special and desirable properties:

1. All of the HL^(k)F(t) components are guaranteed to be “properrotations,” i.e., to have all strictly positive local maxima and allstrictly negative local minima. Again, the need for a proper rotationdecomposition is important in subsequent analysis to reliably determineinstantaneous frequency of the HL^(k)F components.

2. The decomposition is completely determined by the input signal. It isfully adaptive and local in time. Moreover, it is highly computationallyefficient, each step consisting only of comparisons to determineextrema, and followed by a piece-wise linear transformation of buffereddata between successive extrema to produce the desired signalcomponents.

3. The ITD procedure completely preserves temporal information in thesignal. All of the critical points of the HL^(k)F components (i.e., thepoints in time at which their local maxima, local minima, inflectionpoints, etc., occur) coincide precisely with the critical points of theoriginal signal.

4. The ITD method, as exemplarily implemented in the code contained inthe Computer Program Appendix hereof, can be performed in real time tosimultaneously produce output as the information becomes available,thereby eliminating all boundary effects or windowing effects other thaninitial startup effects which last at a particular level only until twoextrema are obtained by the low-frequency component at the prior level,i.e., until full wave information becomes available at these lowerfrequencies. The online implementation allows the method to resolvearbitrarily low frequency information present in the signal and, thus,all frequency information present in the signal.

The ITD method and apparatus of the present invention has broad andimportant applicability to analysis of many different types of signals,including, but not limited to, geophysical signals, seismic data,thermal signals such as sea surface temperature, radiometer signals,environmental signals, biologic signals such as brain waves or heartsignals, proteins, genetic sequences or other data, morphometrics,telecommunications signals, speech or other acoustic signals,crystallographic, spectroscopic, electrical or magnetic signals, objecttrajectories or other physical signals, structural vibrations or othersignals indicative of structural integrity (e.g., movement of bridges orbuildings) such as resonant frequencies of structures, power signalsincluding those in circuits, and signals arising in finance such asstock and bond prices. The method is designed for application to signalsof arbitrary origin, timescale, and/or complexity. It can also be usefulin fusion of data obtained from different sensors, since underlyingcorrelated signals may be uncovered and time-locked together with themethod.

Although the ITD system can incorporate a sifting procedure at eachlevel of decomposition, it is not required. The desirable “properrotation” property of the components, guaranteed by the ITD system, isattained at each level in a single, highly computationally efficient,first step. The only purpose for applying such a sifting procedure wouldbe to smooth the amplitude envelope of the components. While this may bedesirable for certain applications, it must be performed with care, withthe understanding that the sifting process may reduce or even eliminatethe instantaneous information that the ITD otherwise provides withoutany sifting procedure.

Some advantages over the prior art provided by the present inventioninclude:

(a) applicability for analyzing non-stationary and/or non-linearsignals;

(b) ability to operate in real-time due to its computational simplicity,recursive nature, elimination of need for sifting, and absence ofsignificant end effects;

(c) ability to adapt to any timescale and to use complete signalinformation, including all critical points such as inflection points andzero crossings and not just local extrema, thereby allowing weak signalsembedded in stronger signals to be extracted;

(d) ability to provide precise time-frequency-energy localization,either instantaneously via the Hilbert transform or on other timescales,such as the inter-extrema timescale, the inter-critical point timescale,etc.;

(e) ability to interpolate adjacent critical points using the signalitself, sometimes referred to as the self-extractionproperty. The ITDsystem is ‘well-behaved’ even if the data is unevenly sampled or if itis discontinuous. This translates into the ability to preserve intactphase and morphology characteristics of the signal; and

(f) ability to estimate the information content or complexity of thesignal at multiple time and spatial scales.

The characteristics that make the ITD extremely well-suited for theanalysis of nonlinear or non-stationary signals and which fundamentallydifferentiate it from prior art systems are:

(a) The ability to operate in real-time, due to its computationalsimplicity (including the lack of need for sifting, and the uniqueability to generate a proper rotation HF component in one step) and theabsence of significant end effects.

(b) The ability to adapt to any timescale and to use complete signalinformation, including all critical points (e.g., inflection points andzero crossings) and not just local extrema.

(c) The ability to interpolate adjacent critical points using the signalitself (self-extraction property), instead of cubic splines, which causeinstability (overshooting or undershooting) by generating spuriousextrema or by shifting or amplifying existing ones. Additionally, splineinterpolation does not work well when there are sudden changes in thetimescale of the signal, a common phenomenon in analysis ofnon-stationary signals.

(d) The ability to extract single modes in one step, which makes ithighly computationally efficient.

These technical differences translate into important additionaladvantages:

(a) There is absence of significant “end effects.”

(b) The ability to extract weak signals, embedded in stronger onesbecause the ITD scaling is not exclusively given by extrema,

(c) The ITD is “well-behaved” even if the data is unevenly sampled or ifit is discontinuous.

Although the aforementioned advantages, well-suited characteristics, andtechnical differences are described in this iteration sub-section, it isto be understood that these advantages, well-suited characteristics, andtechnical differences are similarly applicable to online ITD, obtainingtime-frequency distribution of signal energy and obtaining instantaneousphase and frequency, analyzing single waves to create new methods offiltering and compression for signals of arbitrary origin, and generaldata analysis as described in the following sub-sections hereof.

B. Online ITD

The ITD method has been recursively implemented for online signalanalysis. Software, written in C and in MATLAB® languages, is providedin the Computer Program Appendix hereof. The online version of thecomputer code is designed to run continuously forward in time and toprocess input signal information as it is obtained, with minimalcomputational expense and minimal time delay. The software detectsextrema in a digital input signal as soon as they occur, and immediatelycomputes the monotonic segment of the corresponding low and highfrequency components on the time interval from the most recent priorextrema to the current, newly found extrema using the procedureexplained herein. The monotonic segment of the LF component is thenconcatenated with the existing signal at the decomposition levelobtained to that point, and the most recent value of the component ischecked to determine whether an extrema has been found at that level. Ifso, the appropriate monotonic segment of the L²F component is thencomputed along with the corresponding segment of the HLF component. Theextrema-triggered process continues down through each level until thelow-frequency monotonic segment simply results in a lengthening of themonotonic segment rather than generating a new extrema at that level.The delay in obtaining information at any given level is simply equal tothe time until the next extrema is obtained at that frequency level.Thus, information is available on the same timescale as that at whichthe input signal is fluctuating, as determined by the time betweensuccessive extrema.

The software requires three signal extrema to begin the decomposition atthe first level. The ability of the method to constrain edge-effects tothe time period prior to the first two extrema allows it to beinterpreted as a start-up transient in information generation, which iscommon in automated signal analysis. From the time of occurrence of thethird extrema in the raw signal forward, or the L^(n)F component beingdecomposed in the case of higher levels of decomposition, the remainingdata at that level is absolutely free from edge-effects.

FIGS. 4A-4D illustrate the procedure followed by the online ITD system.In FIG. 4 a, the input signal to be decomposed is shown in solid lineswith extrema indicated by circles. The baseline component, computed intime up to the local minimum extrema at time t(j−1) is shown in dashedlines. The system is buffering each new input signal value on the solidcurve until such time as a new extrema is detected at time t(j+1). Thedetection of the new extrema triggers computation (via the ITD stepalgorithm) of the monotonic baseline segment for values of t betweent(j−1) and t(j). First, a new baseline node is computed (see solid pointat t(j) in FIG. 4B), and the input signal itself over that interval islinearly transformed to form the monotonic baseline function over thatinterval (see dashed red curve between t(j−1) and t(j) in FIG. 4C). Thelow frequency baseline component segment is then immediately subtractedfrom the input signal to generate the corresponding high frequencycomponent on the [t(j−1), t(j)] interval. The procedure continues, witheach new extrema triggering a similar portion of the decomposition onthe time interval between adjacent extrema. See FIG. 4D for anillustration of resulting baseline and high frequency components forthis example, along with the original signal being decomposed. Notethat: (i) the high frequency component has all positive maxima and allnegative minima (i.e., it is a “proper rotation”), and (ii) the methodis iterated, applying this decomposition procedure to each resultingbaseline component in a manner that is similarly triggered each time thecomponent to be decomposed has a new extrema.

C. Application of the ITD to Obtain the Time-Frequency Distribution ofSignal Energy and New Methods for Obtaining Instantaneous Phase andFrequency.

Application of the ITD to an input signal results in a decompositionthat consists of a set of proper rotation components of successivelylower relative frequencies, along with a monotonic signal trend. Theproper rotations are ideally suited for subsequent analysis to determinethe time-frequency-energy (“TFE”) distribution of the original signal.This can be accomplished, e.g., by using the Hilbert transform if TFEinformation is sought at every data point (i.e., with temporal precisionequal to the signal sampling rate). However, the present invention alsoincludes methods for instantaneous phase angle and instantaneousfrequency computation that do not require use of the Hilbert transform.These methods are “wave-based,” i.e., they define the instantaneousfrequency in a piece-wise manner, each piece corresponding to the timeinterval between successive upcrossing of a proper rotation and usingonly information about the single wave of the proper rotation occurringduring that period and guarantee a monotonically increasing phase anglewhen applied to proper rotation components, whereas numericalcomputations of instantaneous phase using the Hilbert transform do notalways result in this highly desirable property. Instantaneous frequencycan then be obtained from the time derivative of the instantaneous phaseangle. Monotonically increasing phase angles result in instantaneousfrequencies that are never negative. By contrast, in Hilberttransform-based numerical computations of instantaneous phase angle,angle decreases and corresponding negative instantaneous frequenciesoccasionally occur, along with related phase-unwrapping jumps by 2πmultiples. By contrast, the first embodiment of the method describedherein (i.e., the “arcsine approach”) produces instantaneous phaseangles that do not suffer from these problems and, though notnecessarily identical to those obtained via the Hilbert transform,provide a reasonable and useful alternative. The arcsine approachcoincides with the Hilbert transform result for mean-zero trigonometricfunctions and provides more desirable results in cases when the Hilberttransform leads to negative frequencies and/or phase-unwrapping jumps.Moreover, the Hilbert transform approach suffers from boundary effectsthat are overcome by the present invention. One skilled in the art willalso appreciate that while the Hilbert transform and the inventionsdescribed herein may be applied to non-proper rotation signals, theconcept of instantaneous phase angle and instantaneous frequency aremeaningful and avoid ambiguity only when applied to proper rotationsignals.

The first embodiment of the method for instantaneous phase anglecomputation obtains the phase angle, θ(t), from the signal x(t) for onefull wave (i.e., the portion of the signal between successive zeroup-crossings) at a time. The instantaneous phase angle is obtained usingthe arcsine function applied to the positive and negative signalhalf-waves (separated by the zero down-crossing of the wave) afteramplitude normalization, as follows:${\theta(t)} = \left\{ \begin{matrix}{\sin^{- 1}\left( \frac{x(t)}{A_{1}} \right)} & {t \in \left\lbrack {t_{1},t_{2}} \right)} \\{\pi - {\sin^{- 1}\left( \frac{x(t)}{A_{1}} \right)}} & {t \in \left\lbrack {t_{2},t_{3}} \right)} \\{\pi - {\sin^{- 1}\left( \frac{x(t)}{A_{2}} \right)}} & {t \in \left\lbrack {t_{3},t_{4}} \right)} \\{{2\pi} + {\sin^{- 1}\left( \frac{x(t)}{A_{2}} \right)}} & {t \in \left\lbrack {t_{4},t_{5}} \right)}\end{matrix} \right.$where A₁>0 and A₂>0 are the respective amplitudes of the positive andnegative half-waves between the successive zero up-crossings, t₁ and t₅,t₂ is the (first) time of the maxima (A₁) on the positive half-wave, t₃is the zero down-crossing time, and t₄ is the (first) time of the minima(−A₂) on the negative half-wave. Incomplete waves at the beginningand/or end of a proper rotation component, i.e., data prior to the firstzero up-crossing or after the last zero up-crossing, may be treatedaccording to this definition as well over the appropriate sub-intervalof [t₁, t₅]. According to this phase angle definition, the phase angleis zero at every zero up-crossing, π/2 at the local maxima of the wave,π at each zero down-crossing, and 3 π/2 at each local minima. In theon-line ITD, the evolution of the phase angle over any monotonic segmentof a proper rotation can be computed between times of successive extremaas soon as the right-hand extrema is determined and the segment has beenobtained from the ITD decomposition.

The second embodiment of the method for instantaneous phase anglecomputation is similar to the first, in that it is based on theprogression of the proper rotation through each successive full wave,but is designed to provide an even more computationally efficientalternative. In this embodiment, the phase angle is computed for eachwave of the signal as follows: ${\theta(t)} = \left\{ \begin{matrix}{\left( \frac{x(t)}{A_{1}} \right)\frac{\pi}{2}} & {t \in \left\lbrack {t_{1},t_{2}} \right)} \\{{\left( \frac{x(t)}{A_{1}} \right)\frac{\pi}{2}} + {\left( {1 - \frac{x(t)}{A_{1}}} \right)\pi}} & {t \in \left\lbrack {t_{2},t_{3}} \right)} \\{{\left( {- \frac{x(t)}{A_{2}}} \right)\frac{3\pi}{2}} + {\left( {1 + \frac{x(t)}{A_{2}}} \right)\pi}} & {t \in \left\lbrack {t_{3},t_{4}} \right)} \\{{\left( {- \frac{x(t)}{A_{2}}} \right)\frac{3\pi}{2}} + {\left( {1 + \frac{x(t)}{A_{2}}} \right)2\pi}} & {t \in \left\lbrack {t_{4},t_{5}} \right)}\end{matrix} \right.$This approach obtains the instantaneous phase angle via linearinterpolation of the signal value between zero and its value at anextrema. Phase angles that are computed using this second method maydeviate more significantly from those obtained via the Hilbert transformthan those derived from the first embodiment. However, it will beappreciated that either embodiment may be more useful than the other (orHilbert transform-based analysis) for certain applications, such asmeasuring certain types of synchronization between two proper rotations.Which approach is better depends upon such factors as the weights givento physical meaningfulness of the instantaneous frequency and to thecomputational expense of the procedure.

If temporal information on the timescale of the individual waves presentin the decomposed components is all that is required, as opposed toinstantaneous information on the time-scale of the sampling rate asdiscussed above, TFE information may be obtained through quantificationof the size/energy and duration/frequency/wavespeed of the individualwaves of each proper rotation component. Alternatively, one skilled inthe art will appreciate that the ITD method allows the use of othertemporal segmentation for the purpose of TFE distribution construction,such as, e.g., that determined by component critical points, individualmonotonic segments, or component “half-waves” as partitioned byzero-crossings. The results of any of these approaches may be furthersmoothed as desired for visualization and analysis purposes using anynumber of standard techniques well known to those skilled in the art.

FIGS. 5A-5J demonstrate how the ITD may be applied to determine a TFEdistribution of a sample signal and compares its performance to priorart methods of Fourier analysis and wavelet analysis for this task. FIG.5A contains a sample signal. FIGS. 5B-5C contain the TFE distribution ofthe signal obtained via windowed Fourier analysis. Note the rectangulargrid of TFE information naturally derived from this transform and thetrade-offs between good temporal but poor frequency localization (asshown in FIG. 5B), and poor temporal but good frequency localization (asshown in FIG. 5C). Darker regions correspond to those with higher powerthroughout this Figure. It is impossible to simultaneously improve bothtemporal and frequency localization of energy. Quantification isperformed in this Fourier analysis by sequential time windows with apredetermined segmentation of the signal that is not determined in anyway by the signal changes (e.g., local extrema). FIG. 5D illustrates theTFE distribution for the same input signal that is obtained using thefast wavelet transform, sometimes referred to herein as the FWT. Notethe dyadic grid of TFE information naturally derived from thistransform. The FWT has the ability to localize higher frequencyinformation on shorter time scales and lower frequency information onlonger timescales, but still suffers from predetermined temporalsegmentation and inaccuracies in both temporal and frequencylocalization. These are due in part to the bleeding of signal energyacross different levels of resolution in the wavelet transform and areattributed in large part to the temporal and frequency bins that arepredetermined by the choice of basis and not by the signal under study.FIG. 5E shows the ITD-based TFE distribution for the same signal,illustrating the significant advance in time-frequency localizationprovided by the algorithm of the subject invention. Each displayed linesegment's start and end points correspond exactly to the start and endof an actual wave in an ITD proper rotation component and the linesegment shading is determined by the amplitude of the corresponding wave(darker indicates greater amplitude). Since the ITD automaticallydecomposes the input signal into a set of proper rotations and amonotonic trend, the instantaneous phase angles for each proper rotationcomponent can be obtained via the Hilbert transform or the invention'salternative approaches described herein to obtain the correspondinginstantaneous frequencies. FIG. 5F shows the original signal and thefirst four HF components obtained by the ITD. FIG. 5G shows theinstantaneous frequency curves (in light gray when instantaneous poweris insignificant, i.e., below 2×10⁻⁵) corresponding to each of the HFproper rotation components (as shown smoothed and solid black) wheninstantaneous power is non-negligible, i.e., above 0.005, to illustratethis process. The curves in FIG. 5G were obtained using thearcsine-based method for instantaneous frequency computation and weresmoothed using a 0.1 second median filter to enhance visualization. FIG.5H shows the first three proper rotation components obtained byapplication of the ITD to the signal of FIG. 5A. FIGS. 5I and 5J showthe half-wave amplitudes and (smoothed) instantaneous frequencies foreach component obtained via the arcsine-based method for instantaneousphase and instantaneous frequency determination (linestyles are the sameas that of the corresponding proper rotation components). As expected,the instantaneous frequencies of each component decrease at each pointin time as the level increases, but the actual frequency values aredetermined solely by the signal and not predetermined in any manner.

D. Application of the ITD and Single Wave Analysis to Create New Methodsof Filtering and Compression for Signals of Arbitrary Origin.

The ability of the ITD to decompose online any input signal into a setof one or more component signals which have the proper rotation propertymakes it an ideal first step of a two-step process resulting in apowerful and unique new method for nonlinear signal filtering. Inparticular, as each monotonic segment is computed within one or moredecomposition levels, features of the segment can be analyzed andquantified. For example, the amplitude and duration of the segment maybe simply computed from the extrema, i.e., the starting and endingpoints of the segment. One skilled in the art will appreciate that moredetailed features such as the time-weighted feature density of thesegment and associated properties such as its average or median value,its variance, skewness, distance from a template or reference density,etc., may also be quantified for each segment. Individual segments maythen be classified according to the values of theses quantifiedfeatures. A filtered signal can then be constructed by adding together(i.e., superimposing the segments while preserving the time intervals onwhich each occurred) only those segments that satisfy certainconstraints applied to their quantified features. For example, one mayconstruct a signal using only those ITD-created monotonic segments withdurations that, using the data sampling rate, correspond to wavespeedsbetween f_(min) Hz and f_(max) Hz, and with absolute amplitudes of atleast A_(min), e.g., in excess of the 75^(th) percentile of amplitudesfor all segments with the required duration. In this example, the outputcan be essentially interpreted as extracting the waves in the inputsignal that are in the f_(min)−f_(max) frequency band and which havelarger amplitudes than A_(min). Such uses would bear a strongresemblance to the commonly used Fourier-based band-pass filtering, butwould avoid the phase-shifting and waveform distortion drawbacks, whileadding the amplitude differentiation capability. In FIGS. 6A-6G, thisconcept is illustrated when applied to the same input signal used inFIG. 5A (also shown again in FIG. 6A). After the TFE distribution iscomputed, the absolute amplitude and duration/wavespeed of eachindividual wave are obtained. FIGS. 6B-6C show the one-dimensional(marginal) densities of these two features along with verticaldiscriminator in solid lines that were selected for purposes ofillustrating the capabilities of this method. FIG. 6D shows the(absolute amplitude, wavespeed) ordered pairs for each segment.Discriminators of pairs with wavespeed less than 20 Hz, between 20 Hzand 50 Hz, and greater than 50 Hz are applied, along with adiscriminator of amplitude exceeding 0.01 mV. Feature vector pairs(i.e., points in FIG. 6D) falling into certain specific regions of thetwo-dimensional feature range are classified according to the region inwhich they belong. FIG. 6 e displays three filtered componentsconstructed using only those segments with points in specific regions(corresponding to three of the regions in FIG. 6D). Component one, C₁,is the output of a filter that contains the signal reconstructed usingonly those waves that have absolute amplitude exceeding 0.01 mV andwavespeed between 20 Hz and 50 Hz. Component two, C₂, is the output of afilter that contains the signal reconstructed using only those wavesthat have absolute amplitude less than 0.01 mV and wavespeed between 20Hz and 50 Hz. Component three, C₃, is the output of a filter thatcontains the signal reconstructed using only those waves that haveabsolute amplitude exceeding 0.01 mV and wavespeed below 20 Hz. Theoutput obtained by combining some of these components (i.e., C₁+C₂+C₃and C₁+C₂) are shown in FIG. 6F, along with the output of a traditional,linear, [20, 50] Hz. band-pass digital filter provided for comparison.The extrema-preserving properties of this new type of filter are evidentin this FIG. 6F. FIG. 6G again compares an output available by using theITD filter in comparison to the same linear band-pass digital filter,but does so on a slightly larger timescale (10 seconds). The residualsignals obtained by subtracting the filtered outputs from the originalsignal are also shown, to illustrate the ability of the ITD based signalto more completely decompose the raw input signal into what may beinterpreted as abnormal and normal components present in the raw signal.

Since the relative starting time of each segment can also be consideredas a segment feature, temporal relationships between segments can beanalyzed (in addition to other features) to determine which segments toinclude in construction of the filtered signal output. This allows muchmore sophisticated pattern recognition techniques to be incorporatedinto the second step of this filtering process. Of course, there may bemultiple classes of features and an input signal may be decomposed bythis ITD-based nonlinear filter into a multitude of output components.One skilled in the art will also appreciate that adaptation techniquesused in conventional digital filtering may be equally well employed inthis new setting to produce adaptive ITD-based nonlinear filters. Suchtechniques include, for example, cluster-analysis applied to segmentfeatures in order to determine recurring patterns or underlying signaldynamics that may exist in the input signal and which the filter may bedesigned to separate and illuminate. In addition to, or as analternative to single segment-based feature analysis, one may performthe analysis on whole “waves” consisting of pairs of (i.e., consecutive)monotonically increasing and decreasing segments or waves defined byinter-zero-up-crossings of any of the proper rotation components.

Another example that illustrates the power of the invention for signalfiltering and analysis is provided in FIGS. 7A-7H. In these Figures, wedemonstrate how the ITD-based filtering method may be applied todifferentiate between two types of waves that have significantlyoverlapping spectral characteristics and to decompose a raw signal thatis a combination of these two types of waves into components thatpreserve temporal location of extrema and critical points in the rawsignal. Overlap between the power spectral densities of an underlyingsignal and superimposed noise is one of the more difficult problems insignal analysis and detection, with numerous applications. FIG. 7Aillustrates the mix of two alternating signals, one a 9 Hz cosine wavewith an increasing linear trend and the other a 9 Hz sawtooth wave. Thesignals were added to create a raw test signal for use in this example.FIG. 7B shows the power spectral density estimates for each of the twosignals, clearly illustrating the significant overlap in the band around9 Hz. FIG. 7C shows the proper rotation components and monotonic trendobtained via application of ITD to the original combined signal. FIGS.7D-7E illustrate the result of applying single wave analysis to examinevarious features of individual waves obtained by the ITD decomposition.FIG. 7D shows the ratio of individual wave peak to mean signal value,while FIG. 7E shows a plot of wave skewness plotted versus wavekurtosis. A sample discriminator in dashed linestyle is also shown toillustrate the result of feature classification as applied to thesequantified waveform features. FIG. 7F shows the output of the ITD-basedfilter that retains only those waves with peak-to-mean value ratioexceeding 1.65. The corresponding residual signal, i.e., the originalsignal minus the filter output signal, is also shown and illustratesthat peak-to-mean value ratio discrimination of ITD-produced properrotation waves allow reconstruction of a filter output that nearlyperfectly decomposes the original signal into its constituent componentswhile retaining precise temporal information regarding, e.g., extremalocation. FIG. 7G shows the output of the ITD-based filter that retainsonly those waves with skewness-kurtosis pairs above the dasheddiscriminator line, along with the corresponding residual signal. Thisexample illustrates how this new filter method may be easily configuredto decompose an input signal into components with more complex criteria,such as “cutting off the tops of all riding sinusoidal waveforms whileretaining every other signal component” as is done here. FIG. 7H showsthe same signals as in FIG. 7G but zooming in to the time intervalbetween t=9 and t=13 to see details of the filter output.

Another useful capability of the invention involves use of the ITD andsingle wave analysis to provide a means for signal compression. Sincethe ITD proper rotation components preserve extrema and critical pointsof the original signal, single-wave analysis-based compression andsubsequent reconstruction are also able to accurately preserve theextrema and/or critical points deemed significant by the user. The usermay throw away all individual monotonic segments of signal that aredetermined to be insignificant, according to automated analysis of a setof segment features and retain only those segments that are of interestfor later use in reconstruction. Moreover, utilizing cluster analysis ofwaveform shapes, shapes of their derivatives, clustering analysis ofquantifiable signal features, or classification of individual waveformsegments into a small set of ‘cymemes’ or template waves such as, forexample, half sine waves with amplitude and frequency matching thesingle waves, allows for significant compression of non-stationarysignals.

E. Application of the Above Methods to General Data Analysis

One skilled in the art will recognize that the ITD method does notactually require the input data to be a time series or “signal” in thestrict sense of the word, i.e., a sequence of data indexed by time. Themethod is sufficiently flexible to be applied to any set of numericaldata that is a function of a real-valued variable, i.e., one thatsatisfies a “vertical line test.” The index that is typicallyinterpreted herein as a time index could instead just refer to theaforementioned real-valued variable. For example, in any list orsequence of numbers, the variable “t” may be assigned natural numbervalues, i.e., 1, 2, 3, . . . , and interpreted simply as the index intothe list of data. The fact that the ITD does not require uniform spacingof the data in time also allows obvious generalization of the method toother types of data besides time series.

One skilled in the art will also recognize that, just as Fourier-basedfiltering and other decomposition methods can be successfully applied tohigher dimensional analysis, so, too, can the ITD system. This extendsthe above-mentioned application of the ITD to decompose functions ofseveral variables. For example, in image processing applications, onemay be interested in separating an image into a high frequencycomponent, which contains edge transitions between various objects inthe image for example, and a low frequency component, which definebackground colors in the image for example. Two-dimensional waveletand/or Fourier analysis are popular for image processing applications.However, the drawbacks mentioned above for these methods still exist inhigher dimensional analysis. The benefits of the ITD system, such as itsability to preserve extrema and precisely localize time-frequency-energyinformation are preserved in higher dimensional analysis and therebyoffer improvements over the most popular prior art methods. FIGS. 8A-8Cillustrate the application of ITD to decompose a two-dimensional surface(shown in FIG. 8A) into two component surfaces, namely a high frequencysurface (shown in FIG. 8B) and a low frequency surface (shown in FIG.8C). The decomposition in this example is obtained by decomposing eachof the cross-sectional signals of the surface along the grid linesobtained by holding the first independent variable constant, repeatingthe process holding the second variable constant, and averaging the tworesults.

The present invention may include an article of manufacture having acomputer-readable medium comprising code, or a microprocessor havingsuch code embedded therein, wherein such code is capable of causing thecomputer, microprocessor, or other computational apparatus to executethe inventive method.

Further scope of applicability of the present invention will becomeapparent from the detailed description given herein. However, it is tobe understood that the detailed description and specific examples, whileindicating preferred embodiments of the invention, are given by way ofillustration only, since various changes and modifications within thespirit and scope of the invention will become apparent to those skilledin the art from this detailed description. Furthermore, all themathematical expressions are used as a short hand to express theinventive ideas clearly and are not limitative of the claimed invention.

1. (canceled)
 2. (canceled)
 3. (canceled)
 4. (canceled)
 5. (canceled) 6.(canceled)
 7. (canceled)
 8. (canceled)
 9. (canceled)
 10. (canceled) 11.(canceled)
 12. (canceled)
 13. (canceled)
 14. (canceled)
 15. (canceled)16. (canceled)
 17. (canceled)
 18. (canceled)
 19. (canceled) 20.(canceled)
 21. (canceled)
 22. (canceled)
 23. A method of determining atleast one instantaneous phase angle estimate from a proper rotationsignal, including the steps of: (a) inputting a proper rotation signalto a processor; (b) determining a zero-crossing and a local extremum ofthe proper rotation signal; and (c) applying linear interpolation to theproper rotation signal between the zero-crossing and the local extremumto determine an instantaneous phase angle estimate of the properrotation signal.
 24. A method of determining at least one instantaneousphase angle estimate from a proper rotation signal, including the stepsof: (a) inputting a proper rotation signal to a processor; (b)extracting an amplitude-normalized half wave from the proper rotationsignal; and (c) applying an arcsine function to the amplitude-normalizedhalf wave to determine an instantaneous phase angle estimate of theproper rotation signal.
 25. A method of determining at least oneinstantaneous amplitude estimate from a proper rotation signal,including the steps of: (a) inputting a proper rotation signal to aprocessor; (b) determining a zero-crossing and a local extremum of theproper rotation signal; and (c) determining an instantaneous amplitudeestimate on the interval between the zero-crossing and the localextremum, using the value of the proper rotation at the local extremum.26. A method of determining time-frequency-energy information from asignal, comprising the steps of: (a) decomposing the signal into a setof one or more proper rotation signals and a residual baseline trendsignal; (b) determining at least one zero-crossing and at least onelocal extrema for the at least one proper rotation signals; (c)determining at least one time-frequency-energy information feature ofthe at least one proper rotation signal, on the interval between the atleast one zero-crossing and the at least one local extrema, wherein theat least one feature is selected from the group comprising aninstantaneous amplitude estimate, an instantaneous energy estimate, aninstantaneous phase angle estimate, and an instantaneous frequencyestimate.